##########julia file for creating simulated data from model
function Simulate(prim::Primitives, est::Estimands, res::Results, simul_data::Array{Float64,2}, lhoods::Array{Float64, 2}; race::Int64 = 1)
    @unpack μ_grid, μs_grid, e_grid, m_grid, p_grid, ℓ_grid, a_grid_1, x_grid_1, ac_grid, f_grid, ℓp_type_grid, h_grid, e_s_grid = prim #grids
    @unpack n_μ, n_μs, n_e, n_m, n_p, n_ℓ, n_a_1, n_x_1, n_ac, n_f, n_ℓp, n_h, n_e_s = prim #grid dimensions
    @unpack β, λ_0s_vec, λ_1s_vec, λ_2s_vec, λ_3s_vec, div_chars = prim #discount factor and spouse wage parameters
    @unpack α_1, α_2, α_3, α_4, α_5, α_6, α_7, α_8, α_9, α_10, τ_s, τ_p, τ_p_2, λ_0, λ_1, λ_2, λ_3, λ_4, λ_5, γ_0, γ_1, γ_2, γ_3, γ_4, σ_ε, σ_ξ = est
    @unpack γ_5, α_μ, α_x, τ_prob, α_amen_1, α_amen_2, α_amen_3, λ_6, μ_prob = est
    @unpack emax_1, emax_2, cutoff_1, cutoff_1_fert = res
    @unpack θ_1, θ_2, θ_3, θ_4, σ_θ, θ_5 = est

    #yank out appropriate spousal wage process
    λ_0s = λ_0s_vec[race]
    λ_1s = λ_1s_vec[race]
    λ_2s = λ_2s_vec[race]
    λ_3s = λ_3s_vec[race]

    dist_fert = Normal(0, σ_θ)
    nrow = length(simul_data[:,1]) #store number of rows in data
    Random.seed!(1234)
    nrows_simul, ncol = 0, 21

    #figure out how many rows we want
    for s = 1:10
        for r = 1:nrow
            row = simul_data[r,:] #access particular row in data; also a preallocation
            i_a_start = Int64(row[8])
            i_a_max = Int64(row[16]) #highest age individual is observed in data, to deal with non-random attrition in sample
            for i = i_a_start:i_a_max
                nrows_simul += 1
            end
        end
    end
    output = zeros(nrows_simul, ncol) #vector to fill

    #index of simulated row we're on
    i_row = 0
    for s = 1:10
        for r = 1:nrow #loop over individuals
            row = simul_data[r,:] #access particular row in data; also a preallocation

            #tyep probabilities
            p_τ = 1-τ_prob #probbility of type 1s types
            p_μ = 1-μ_prob

            #add type probabilities so sthat we're doing Bayes rule correctly
            lhood_row = copy(lhoods[r,:]) #make sure we don't mess up with mutability stuff
            #lhood_row = ones(6) #make sure we don't mess up with mutability stuff
            lhood_row[1] *= (p_τ * p_μ)
            lhood_row[2] *= (p_τ * (1-p_μ))
            lhood_row[3] *= ((1-p_τ) * p_μ)
            lhood_row[4] *= ((1-p_τ) * (1-p_μ))
            

            #draw types
            type_probs = lhood_row./sum(lhood_row)
            type_draw = rand() #uniform 0-1 draw
            type = 1

            for i = 1:3
                if type_draw>sum(type_probs[1:i]) && type_draw<(sum(type_probs[1:i+1]))
                    type = i + 1
                end
            end

            #determine type components
            i_τ = 1
            if type>=3
                i_τ = 2
            end
            
            i_μ = 1
            if type == 2 || type == 4
                i_μ = 2
            end

            #not important stuff
            id, year = Int64(row[1]), Int64(row[2])

            #obtain starting indices
            i_e, i_m, i_p, i_ℓ, i_a, i_x, i_ac, i_ℓp, i_e_s =
            Int64(row[4]), Int64(row[5]), Int64(row[6]), Int64(row[7]), Int64(row[8]), Int64(row[9]),
            Int64(row[10]), Int64(row[12]), Int64(row[19])
            i_a_start = Int64(row[8])
            i_a_max = Int64(row[16]) #highest age individual is observed in data, to deal with non-random attrition in sample

            if i_e_s == 0
                i_e_s = 1
            end

            #begin main simulation loop
            for i = i_a_start:i_a_max #simulate until final age in data
                i_row+=1
                #unpack some useful stuff
                e, age, ac, m, x, μ = e_grid[i_e], i_a + 21, ac_grid[i_ac], m_grid[i_m], x_grid_1[i_x], μ_grid[i_μ]
                e_s = e_s_grid[i_e_s]

                #nab location wage fixed effect
                η_h, δ, coli, η = div_chars[i_ℓp, 2], div_chars[i_ℓp, 5], div_chars[i_ℓp, 4], div_chars[i_ℓp, 9]
                d2shore, amen_index, cooling = div_chars[i_ℓp, 6], div_chars[i_ℓp, 7], div_chars[i_ℓp, 8]

                if i_ℓ!=1
                    η_h, δ, coli, η = div_chars[i_ℓ-1, 2], div_chars[i_ℓ-1, 5], div_chars[i_ℓ-1, 4], div_chars[i_ℓ-1, 9]
                    d2shore, amen_index, cooling = div_chars[i_ℓ-1, 6], div_chars[i_ℓ-1, 7], div_chars[i_ℓ-1, 8]
                end

                Δ = γ_0 - γ_1 * e + γ_2 * (age-21) + γ_3 * (ac>-1) + γ_4 * (m>0) #compute moving cost
                wage_predict = λ_0 + λ_1*e + λ_2*x + λ_3*x^2 + μ + η * (1 + λ_6*e) + λ_4 * (ac==0 || ac==1) + λ_5 * (ac>1)
                μs = μs_grid[i_m]
                xh = age-(18 + 4*e_s) #experience of husband
                w_s = exp(λ_0s + λ_1s*e_s + λ_2s*xh + λ_3s*xh^2 + μs + η_h) * (m>0) #spouse wsages: 0 if unmarried

                #determine fertility decision
                cutoff_fert = cutoff_1_fert[i_μ, i_e, i_m, i_p, i_ℓ, i_a, i_x, i_ac, i_ℓp, i_τ, i_e_s]
                draw_f = rand(dist_fert)
                f, i_f = 0, 1
                util_fert = 0 #utiltiy associated with fertiltiy decision
                if draw_f>cutoff_fert
                    f, i_f = 1, 2 
                    util_fert = (θ_1 - θ_2 * (m>0) + θ_3 * age + θ_4 * age * (m>0) - θ_5 * e - draw_f) * -1
                end

                trans_states = [e, m, age, e_s] #things that govern probabilities of marriage, fertility
                i_ac_p = 1 #index of next-period child age state
                if ac>=0 && ac<4 #kid aged  0 to 3
                    i_ac_p = i_ac+1 #next period: kid aged one year older
                end

                if f == 1 #check fertility status
                    i_ac_p = 2 #kdi aged 0 next period if woman is pregnant
                end

                #determine hours decision
                cutoff = cutoff_1[i_μ, i_e, i_m, i_p, i_ℓ, i_a, i_x, i_ac, i_f, i_ℓp, i_τ, i_e_s]
                draw = rand(Normal(0, σ_ε))
                draw_ξ = rand(Normal(0, σ_ξ))

                wage, h = 0, 0

                if draw>cutoff
                    wage = exp(wage_predict + draw + draw_ξ) #add in measuremet error
                    h = 1
                end

                #utility, part 1
                α_1p, α_2p, α_3p = α_1, α_2, α_3
                util_bonus = α_4 * (i_ℓ == 1)
                if i_ac>1
                    α_1p, α_2p, α_3p = α_5, α_6, α_7
                    util_bonus = α_8 * (i_ℓ == 1)
                end
                α_9p = α_9/coli
                α_1p /= coli

                #transfers
                τ_p_star = τ_p
                if m>0
                    τ_p_star = τ_p_2
                end

                if i_τ == 2
                    τ_p_star = 0.0
                end

                ccc = δ * max(0, 1 - τ_s*(m>0) - τ_p_star * (i_ℓ == 1)) * (ac>-1) #childcare costs
                c = wage + w_s - ccc
                util = α_1p * c + (1-h) * (α_2p + α_μ * (i_μ-1) + α_x * x + α_9p * w_s + α_10 * e) + α_3p * (h+1!=i_p)
                util += (util_bonus + α_amen_1 * d2shore + α_amen_2 * amen_index + α_amen_3 * cooling)
                util += util_fert

                i_p_p = h + 1 #next-period past-participation state
                i_x_p = min(i_x + h, n_x_1) #next-period experience
                i_x_p_2 = i_x + h

                #get moving probabilities
                emax_grid = zeros(n_ℓ+1)
                if i_a == n_a_1 #in terminal age of phase 1
                    #parent location option
                    emax_grid[1] = β * emax_2[i_μ, i_e, i_m, i_p_p, 1, 1, i_x_p_2, i_ac_p, i_ℓp, i_τ, i_e_s] - (Δ + γ_5 * div_chars[i_ℓp, 3]) * (i_ℓ!=1)

                    #moving options
                    for i_ℓc = 2:(n_ℓ)
                        emax_grid[i_ℓc] = β * emax_2[i_μ, i_e, i_m, i_p_p, i_ℓc, 1, i_x_p_2, i_ac_p, i_ℓp, i_τ, i_e_s] - (Δ + γ_5 * div_chars[i_ℓc-1, 3])
                    end

                    #stay8ing option
                    emax_grid[n_ℓ+1] = β * emax_2[i_μ, i_e, i_m, i_p_p, i_ℓ, 1, i_x_p_2, i_ac_p, i_ℓp, i_τ, i_e_s]
                elseif i_a!=n_a_1 #not in terminal age

                      #loop over potential marriage/fertility realizations
                      for i_mp = 1:n_m, i_e_s_p = 1:n_e_s
                        mp = m_grid[i_mp] #standardize
                        e_s_p = e_s_grid[i_e_s_p]
                        prob = Compute_Prob_Marr(mp, e_s_p, trans_states, prim) #compute probability of realization. Slight bottleneck here.
                        emax_grid[1] += (β * emax_1[i_μ, i_e, i_mp, i_p_p, 1, i_a+1, i_x_p, i_ac_p, i_ℓp, i_τ, i_e_s_p] - (Δ + γ_5 * div_chars[i_ℓp, 3]) * (i_ℓ!=1)) * prob
                        for i_ℓc = 2:(n_ℓ)
                            emax_grid[i_ℓc] += (β * emax_1[i_μ, i_e, i_mp, i_p_p, i_ℓc, i_a+1, i_x_p, i_ac_p, i_ℓp, i_τ, i_e_s_p] - (Δ + γ_5 * div_chars[i_ℓc-1, 3])) * prob
                        end
                        emax_grid[n_ℓ+1] += (β * emax_1[i_μ, i_e, i_mp, i_p_p, i_ℓ, i_a+1, i_x_p, i_ac_p, i_ℓp, i_τ, i_e_s_p]) * prob
                    end
                end

                ℓ_shocks = rand(Gumbel(0, 1), n_ℓ+1)
                lp_choice = maximum(emax_grid .+ ℓ_shocks)
                lp = findmin(abs.(emax_grid.+ℓ_shocks.-lp_choice))[2] #index of choice
                if lp == 1 && i_ℓ == 1 #opting to stay in parent location
                    lp = 11 #recode to stay
                end
                if i == i_a_max
                    lp = 99 #location nob observed in final row of data
                end

                #utility, part 2
                if lp!=99
                    util += ℓ_shocks[lp] #utiltiy syhock associated with chosen location
                end

                if lp == 1 #moved to parent location
                    util -= (Δ + γ_5 * div_chars[i_ℓp, 3])
                elseif lp>2 && lp<11 #moved somewhere
                    util -= (Δ + γ_5 * div_chars[lp-1, 3])
                end

                #write line of data
                line = [id year i_μ i_e i_m i_p i_ℓ i_a i_x i_ac i_f i_ℓp h wage lp s i_a_max i_τ util c/coli e_s]
                output[i_row,:] = line
                year += 1

                ####determine next-period indices
                #marriage/fertility
                draw_m = rand()
                prob_m = [Compute_Prob_Marr(0, 0, trans_states, prim), Compute_Prob_Marr(1, 0, trans_states, prim), Compute_Prob_Marr(2, 0, trans_states, prim), 
                Compute_Prob_Marr(0, 1, trans_states, prim), Compute_Prob_Marr(1, 1, trans_states, prim), Compute_Prob_Marr(2, 1, trans_states, prim)]
                i_m, i_e_s = 3, 2 #preallocate
                
                if draw_m<prob_m[1]
                    i_m, i_e_s = 1, 1
                elseif draw_m>prob_m[1] && draw_m<sum(prob_m[1:2])
                    i_m, i_e_s = 2, 1
                elseif draw_m>sum(prob_m[1:2]) && draw_m<sum(prob_m[1:3])
                    i_m, i_e_s = 3, 1
                elseif draw_m>sum(prob_m[1:3]) && draw_m<sum(prob_m[1:4])
                    i_m, i_e_s = 1, 2
                elseif draw_m>sum(prob_m[1:4]) && draw_m<sum(prob_m[1:5])
                    i_m, i_e_s = 2, 2
                end

                #past participation, experience, age, age of youngest child
                i_p = i_p_p
                i_x = i_x_p
                i_a+=1
                i_ac = i_ac_p

                #location
                if lp!=11 #moving
                    i_ℓ = lp #update
                end
            end
        end
    end
    output #return simulated data
end


##########function to simulate entire lifecycle for individuals
function Simulate_lifecycle(prim::Primitives, est::Estimands, res::Results, simul_data::Array{Float64,2}, lhoods::Array{Float64, 2}; cfact::Int64 = 0, race::Int64 = 0)
    @unpack μ_grid, μs_grid, e_grid, m_grid, p_grid, ℓ_grid, a_grid_1, x_grid_1, ac_grid, f_grid, ℓp_type_grid, h_grid, e_s_grid = prim #grids
    @unpack x_grid_2, a_grid_2, x_grid_3, a_grid_3, n_e_s = prim
    @unpack n_μ, n_μs, n_e, n_m, n_p, n_ℓ, n_a_1, n_x_1, n_ac, n_f, n_ℓp, n_h, n_a_2, n_a_3, n_x_2, n_x_3 = prim #grid dimensions
    @unpack β, λ_0s_vec, λ_1s_vec, λ_2s_vec, λ_3s_vec, div_chars = prim #discount factor and spouse wage parameters
    @unpack α_1, α_2, α_3, α_4, α_5, α_6, α_7, α_8, α_9, α_10, τ_s, τ_p, τ_p_2, λ_0, λ_1, λ_2, λ_3, λ_4, λ_5, γ_0, γ_1, γ_2, γ_3, γ_4, σ_ε, σ_ξ = est
    @unpack γ_5, α_μ, α_x, τ_prob, α_amen_1, α_amen_2, α_amen_3, λ_6, μ_prob  = est
    @unpack emax_1, emax_2, emax_3, cutoff_1, cutoff_2, cutoff_3, cutoff_1_fert = res
    @unpack θ_1, θ_2, θ_3, θ_4, σ_θ, θ_5 = est

    λ_0s = λ_0s_vec[race]
    λ_1s = λ_1s_vec[race]
    λ_2s = λ_2s_vec[race]
    λ_3s = λ_3s_vec[race]

    dist_fert = Normal(0, σ_θ)
    nrow = length(simul_data[:,1]) #store number of rows in data
    Random.seed!(1234)

    nrows_simul, ncol = 0, 24

    #figure out how many rows we want
    for s = 1:10
        for r = 1:nrow
            row = simul_data[r,:] #access particular row in data; also a preallocation
            i_a_start = Int64(row[8])
            for i = i_a_start:44
                nrows_simul += 1 #22 through 65
            end
        end
    end
    output = zeros(nrows_simul, ncol) #vector to fill

    #index of simulated row we're on
    i_row = 0
    for s = 1:10
        for r = 1:nrow #loop over individuals
            row = simul_data[r,:] #access particular row in data; also a preallocation

            #tyep probabilities
            p_τ = 1-τ_prob #probbility of type 1s types
            p_μ = 1-μ_prob

            #add type probabilities so sthat we're doing Bayes rule correctly
            lhood_row = copy(lhoods[r,:]) #make sure we don't mess up with mutability stuff
            #lhood_row = ones(6) #make sure we don't mess up with mutability stuff
            lhood_row[1] *= (p_τ * p_μ)
            lhood_row[2] *= (p_τ * (1-p_μ))
            lhood_row[3] *= ((1-p_τ) * p_μ)
            lhood_row[4] *= ((1-p_τ) * (1-p_μ))
            

            #draw types
            type_probs = lhood_row./sum(lhood_row)
            type_draw = rand() #uniform 0-1 draw
            type = 1

            for i = 1:3
                if type_draw>sum(type_probs[1:i]) && type_draw<(sum(type_probs[1:i+1]))
                    type = i + 1
                end
            end

            #determine type components
            i_τ = 1
            if type>=3
                i_τ = 2
            end
            
            i_μ = 1
            if type == 2 || type == 4
                i_μ = 2
            end
          
            #not important stuff
            id, year = Int64(row[1]), Int64(row[2])

            #obtain starting indices
            i_e, i_m, i_p, i_ℓ, i_a, i_x, i_ac, i_ℓp, i_e_s =
            Int64(row[4]), Int64(row[5]), Int64(row[6]), Int64(row[7]), Int64(row[8]), Int64(row[9]),
            Int64(row[10]), Int64(row[12]), Int64(row[19])
            i_a_start = Int64(row[8])
            phase = 1 #phase of lifecycle
            #i_μ = Int64(row[3]) #if μ obsereved, otherwise leave commented out

            if i_e_s == 0
                i_e_s = 1
            end

            #original η
            η_orig = div_chars[i_ℓp, 9]

            if i_ℓ!=1
                η_orig = div_chars[i_ℓ-1, 9]
            end

            util_start = emax_1[i_μ, i_e, i_m, i_p, i_ℓ, 1, i_x, i_ac, i_ℓp, i_τ, i_e_s]

            #begin main simulation loop
            for i = i_a_start:44 #simulate until final age in data
                i_row+=1

                #unpack some useful stuff
                e, age, ac, m, μ, = e_grid[i_e], i_a + 21, ac_grid[i_ac], m_grid[i_m], μ_grid[i_μ]
                e_s = e_s_grid[i_e_s]


                x = 0
                if phase == 1
                    x = x_grid_1[i_x]
                elseif phase == 2
                    age, x = i_a + 40, x_grid_2[i_x]
                elseif phase == 3
                    age, x = i_a + 44, x_grid_3[i_x]
                end

                #nab location wage fixed effect
                η_h, δ, coli, η = div_chars[i_ℓp, 2], div_chars[i_ℓp, 5], div_chars[i_ℓp, 4], div_chars[i_ℓp, 9]
                d2shore, amen_index, cooling = div_chars[i_ℓp, 6], div_chars[i_ℓp, 7], div_chars[i_ℓp, 8]

                if i_ℓ!=1
                    η_h, δ, coli, η = div_chars[i_ℓ-1, 2], div_chars[i_ℓ-1, 5], div_chars[i_ℓ-1, 4], div_chars[i_ℓ-1, 9]
                    d2shore, amen_index, cooling = div_chars[i_ℓ-1, 6], div_chars[i_ℓ-1, 7], div_chars[i_ℓ-1, 8]
                end

                ####check for counterfactuals                
                if cfact == 2
                    δ = 0 #national subsidy
                elseif cfact == 3 && i_ℓ == 1 
                    δ = 0 #local subsidy
                elseif cfact == 4
                    δ/=2
                elseif cfact == 5
                    δ -= 0.25
                end

                if i_m>1 && e == 1
                    δ*=1.1
                else
                    δ*=0.9
                end

                Δ = γ_0 - γ_1 * e + γ_2 * (age-21) + γ_3 * (ac>-1) + γ_4 * (m>0)    #compute moving cost
                wage_predict = λ_0 + λ_1*e + λ_2*x + λ_3*x^2 + μ + η * (1 + λ_6*e) + λ_4 * (ac==0 || ac==1) + λ_5 * (ac>1)
                μs = μs_grid[i_m]
                xh = age-(18 + 4*e_s) #experience of husband
                w_s = exp(λ_0s + λ_1s*e_s + λ_2s*xh + λ_3s*xh^2 + μs + η_h) * (m>0) #spouse wsages: 0 if unmarried
                
                #determine fertility decision
                i_f, f = 1, 0
                util_fert = 0 #utiltiy associated with fertiltiy decision

                if age < 41 #assuming after 35, no fertility
                    #determine fertility decision
                    cutoff_fert = cutoff_1_fert[i_μ, i_e, i_m, i_p, i_ℓ, i_a, i_x, i_ac, i_ℓp, i_τ, i_e_s]
                    draw_f = rand(dist_fert)

                    if draw_f>cutoff_fert
                        f, i_f = 1, 2 
                        util_fert = (θ_1 - θ_2 * (m>0) + θ_3 * age + θ_4 * age * (m>0) - θ_5 * e - draw_f) * -1
                    end

                end
                
                trans_states = [e, m, age, e_s] #things that govern probabilities of marriage, fertility
                i_ac_p = 1 #index of next-period child age state
                if ac>=0 && ac<4 #kid aged  0 to 3
                    i_ac_p = i_ac+1 #next period: kid aged one year older
                end
                if f == 1 #check fertility status
                    i_ac_p = 2 #kdi aged 0 next period if woman is pregnant
                end

                #determine hours decision
                cutoff = 0
                if phase == 1
                    cutoff = cutoff_1[i_μ, i_e, i_m, i_p, i_ℓ, i_a, i_x, i_ac, i_f, i_ℓp, i_τ, i_e_s]
                elseif phase == 2
                    cutoff = cutoff_2[i_μ, i_e, i_m, i_p, i_ℓ, i_a, i_x, i_ac, i_ℓp, i_τ, i_e_s]
                elseif phase == 3
                    cutoff = cutoff_3[i_μ, i_e, i_m, i_p, i_ℓ, i_a, i_x, i_ℓp, i_e_s]
                end

                draw = rand(Normal(0, σ_ε))
                draw_ξ = rand(Normal(0, σ_ξ))
                wage, h = 0, 0
                wage_nomove = 0

                if draw>cutoff
                    wage = exp(wage_predict + draw + draw_ξ) #add in measuremet error
                    wage_nomove = exp(wage_predict + draw + draw_ξ - η + η_orig) #add in measuremet error
                    h = 1
                end

                #utility, part 1
                α_1p, α_2p, α_3p = α_1, α_2, α_3
                util_bonus = α_4 * (i_ℓ == 1)
                if i_ac>1
                    α_1p, α_2p, α_3p = α_5, α_6, α_7
                    util_bonus = α_8 * (i_ℓ == 1)
                end
                α_9p = α_9/coli
                α_1p /= coli

                #transfers
                τ_p_star = τ_p
                if m>0
                    τ_p_star = τ_p_2
                end

                if i_τ == 2
                    τ_p_star = 0.0
                end

                ccc = δ * max(0, 1 - τ_s*(m>0) - τ_p_star * (i_ℓ == 1)) * (ac>-1) #childcare costs

                if cfact == 1 
                    ccc = δ * max(0, 1 - τ_s*(m>0) - τ_p_star) * (ac>-1) #grandparent transfer happens regardless of location
                end

                c = wage + w_s - (ccc * h)
                util = α_1p * c + (1-h) * (α_2p + α_μ * (i_μ-1) + α_x * x + α_9p * w_s + α_10 * e) + α_3p * (h+1!=i_p)
                util += (util_bonus + α_amen_1 * d2shore + α_amen_2 * amen_index + α_amen_3 * cooling)
                util += util_fert

                i_p_p = h + 1 #next-period past-participation state
                i_x_p = min(i_x + h, n_x_1) #next-period experience
                if phase == 2
                    i_x_p = min(i_x + h, n_x_2) #next-period experience
                elseif phase == 3
                    i_x_p = min(i_x + h, n_x_3) #next-period experience
                end
                i_x_p_2 = i_x + h #if jumping a phase next period

                #get moving probabilities
                emax_grid = zeros(n_ℓ+1)
                if phase == 1
                    if i_a == n_a_1 #in terminal age of phase 1
                        #parent location option
                        emax_grid[1] = β * emax_2[i_μ, i_e, i_m, i_p_p, 1, 1, i_x_p_2, i_ac_p, i_ℓp, i_τ, i_e_s] -  (Δ + γ_5 * div_chars[i_ℓp, 3]) * (i_ℓ!=1)

                        #moving options
                        for i_ℓc = 2:(n_ℓ)
                            emax_grid[i_ℓc] = β * emax_2[i_μ, i_e, i_m, i_p_p, i_ℓc, 1, i_x_p_2, i_ac_p, i_ℓp, i_τ, i_e_s] -  (Δ + γ_5 * div_chars[i_ℓc-1, 3])
                        end

                        #stay8ing option
                        emax_grid[n_ℓ+1] = β * emax_2[i_μ, i_e, i_m, i_p_p, i_ℓ, 1, i_x_p_2, i_ac_p, i_ℓp, i_τ, i_e_s]
                    elseif i_a!=n_a_1 #not in terminal age
                        
                        
                        #loop over potential marriage/fertility realizations
                        for i_mp = 1:n_m, i_e_s_p = 1:n_e_s
                            mp = m_grid[i_mp] #standardize
                            e_s_p = e_s_grid[i_e_s_p]
                            prob = Compute_Prob_Marr(mp, e_s_p, trans_states, prim) #compute probability of realization. Slight bottleneck here.
                            emax_grid[1] += (β * emax_1[i_μ, i_e, i_mp, i_p_p, 1, i_a+1, i_x_p, i_ac_p, i_ℓp, i_τ, i_e_s_p] - (Δ + γ_5 * div_chars[i_ℓp, 3]) * (i_ℓ!=1)) * prob
                            for i_ℓc = 2:(n_ℓ)
                                emax_grid[i_ℓc] += (β * emax_1[i_μ, i_e, i_mp, i_p_p, i_ℓc, i_a+1, i_x_p, i_ac_p, i_ℓp, i_τ, i_e_s_p] - (Δ + γ_5 * div_chars[i_ℓc-1, 3])) * prob
                            end
                            emax_grid[n_ℓ+1] += (β * emax_1[i_μ, i_e, i_mp, i_p_p, i_ℓ, i_a+1, i_x_p, i_ac_p, i_ℓp, i_τ, i_e_s_p]) * prob
                        end
                    end
                elseif phase == 2
                    if i_a == n_a_2 #in terminal age of phase 1
                        #parent location option
                        emax_grid[1] = β * emax_3[i_μ, i_e, i_m, i_p_p, 1, 1, i_x_p_2, i_ℓp, i_e_s] -  (Δ + γ_5 * div_chars[i_ℓp, 3]) * (i_ℓ!=1)

                        #moving options
                        for i_ℓc = 2:(n_ℓ)
                            emax_grid[i_ℓc] = β * emax_3[i_μ, i_e, i_m, i_p_p, i_ℓc, 1, i_x_p_2, i_ℓp, i_e_s] -  (Δ + γ_5 * div_chars[i_ℓc-1, 3])
                        end

                        #stay8ing option
                        emax_grid[n_ℓ+1] = β * emax_3[i_μ, i_e, i_m, i_p_p, i_ℓ, 1, i_x_p_2, i_ℓp, i_e_s]

                    elseif i_a!=n_a_2 #not in terminal age

                        emax_grid[1] = β * emax_2[i_μ, i_e, i_m, i_p_p, 1, i_a + 1, i_x_p_2, i_ac_p, i_ℓp, i_τ, i_e_s] -  (Δ + γ_5 * div_chars[i_ℓp, 3]) * (i_ℓ!=1)

                        #moving options
                        for i_ℓc = 2:(n_ℓ)
                            emax_grid[i_ℓc] = β * emax_2[i_μ, i_e, i_m, i_p_p, i_ℓc, i_a + 1, i_x_p_2, i_ac_p, i_ℓp, i_τ, i_e_s] -  (Δ + γ_5 * div_chars[i_ℓc-1, 3])
                        end

                        #stay8ing option
                        emax_grid[n_ℓ+1] = β * emax_2[i_μ, i_e, i_m, i_p_p, i_ℓ, i_a + 1, i_x_p_2, i_ac_p, i_ℓp, i_τ, i_e_s]
                    end
                elseif phase == 3
                    if i_a!=n_a_3 #not in terminal age

                        emax_grid[1] = β * emax_3[i_μ, i_e, i_m, i_p_p, 1, i_a + 1, i_x_p, i_ℓp, i_e_s] -  (Δ + γ_5 * div_chars[i_ℓp, 3]) * (i_ℓ!=1)

                        #moving options
                        for i_ℓc = 2:(n_ℓ)
                            emax_grid[i_ℓc] = β * emax_3[i_μ, i_e, i_m, i_p_p, i_ℓc, i_a + 1, i_x_p, i_ℓp, i_e_s] -  (Δ + γ_5 * div_chars[i_ℓc-1, 3])
                        end

                        #stay8ing option
                        emax_grid[n_ℓ+1] = β * emax_3[i_μ, i_e, i_m, i_p_p, i_ℓ, i_a + 1, i_x_p, i_ℓp, i_e_s]
                    end
                end

                ℓ_shocks = rand(Gumbel(0, 1), n_ℓ+1)
                lp_choice = maximum(emax_grid .+ ℓ_shocks)
                lp = findmin(abs.(emax_grid.+ℓ_shocks.-lp_choice))[2] #index of choice
                if lp == 1 && i_ℓ == 1 #opting to stay in parent location
                    lp = 11 #recode to stay
                end
                
                if i == 44
                    lp = 99 #location nob observed in final row of data
                end

                #utility, part 2
                if lp!=99
                    util += ℓ_shocks[lp] #utiltiy syhock associated with chosen location
                end

                if lp == 1 #moved to parent location
                    util -= (Δ + γ_5 * div_chars[i_ℓp, 3])
                elseif lp>2 && lp<11 #moved somewhere
                    util -= (Δ + γ_5 * div_chars[lp-1, 3])
                end

                #write line of data
                line = [id year i_μ i_e i_m i_p i_ℓ i_a i_x i_ac i_f i_ℓp h wage lp s age phase Δ mean(emax_grid) i_τ util util_start e_s]
                output[i_row,:] = line
                year += 1

                ####determine next-period indices
                #marriage/fertility
                draw_m = rand()
                prob_m = [Compute_Prob_Marr(0, 0, trans_states, prim), Compute_Prob_Marr(1, 0, trans_states, prim), Compute_Prob_Marr(2, 0, trans_states, prim), 
                Compute_Prob_Marr(0, 1, trans_states, prim), Compute_Prob_Marr(1, 1, trans_states, prim), Compute_Prob_Marr(2, 1, trans_states, prim)]
                i_m, i_e_s = 3, 2 #preallocate
                
                if draw_m<prob_m[1]
                    i_m, i_e_s = 1, 1
                elseif draw_m>prob_m[1] && draw_m<sum(prob_m[1:2])
                    i_m, i_e_s = 2, 1
                elseif draw_m>sum(prob_m[1:2]) && draw_m<sum(prob_m[1:3])
                    i_m, i_e_s = 3, 1
                elseif draw_m>sum(prob_m[1:3]) && draw_m<sum(prob_m[1:4])
                    i_m, i_e_s = 1, 2
                elseif draw_m>sum(prob_m[1:4]) && draw_m<sum(prob_m[1:5])
                    i_m, i_e_s = 2, 2
                end

                #past participation, experience, age, age of youngest child
                i_p = i_p_p
                i_x = i_x_p
                i_a+=1
                i_ac = i_ac_p

                #check to see if we're moving to next phase next period
                if age==40
                    phase = 2
                    i_a = 1 #reset age index
                elseif age==44
                    phase = 3
                    i_a = 1 #reset age index
                end

                #location
                if lp!=11 #moving
                    i_ℓ = lp #update
                end
            end
        end
    end
    output #return simulated data
end














############
